/*	
	This file is part of the Snoopy code.

    Snoopy code is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Snoopy code is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Snoopy code.  If not, see <http://www.gnu.org/licenses/>.
*/



#ifndef __SNOOPY_H__
#define __SNOOPY_H__

#include <math.h>
#include <complex.h>
#include <fftw3.h>

#ifdef MPI_SUPPORT
#include <mpi.h>
#include <fftw3-mpi.h>
#endif

#ifdef _OPENMP
#include <omp.h>
#endif

#include "gvars.h"
#include "error.h"

// Enforce No elssaser when not using MHD
#ifdef ELSASSER_FORMULATION
#ifndef MHD
#undef ELSASSER_FORMULATION
#endif
#endif

#ifdef MPI_SUPPORT
#define		MPI_Printf			if (g.main) printf
#else
#define		MPI_Printf			printf
#endif

// Filenames

#define	OUTPUT_DUMP					"dump.dmp"			/**< Dump files filename. */
#define OUTPUT_DUMP_SAV				"dump_sav.dmp"      /**< Previous (saved) output dump. */
#define OUTPUT_DUMP_WRITE			"dump_write.dmp"	/**< dump currently written. */

#define		IDX3D			(k + j * g.nzloc + g.nzloc * g.nyloc * i)  /**< General wrapper for 3D Arrays */

// Structures

#ifdef WITH_PARTICLES
struct Particle	{
	double x;
	double y;
	double z;
	double vx;
	double vy;
	double vz;
	double mass;	// unused for the moment
	double stime;	// Stopping time
};
#endif

struct Field {
	int	   nfield;
	double complex **farray;
	char   **fname;
	
	 // This are the actual fields. 
	double complex *vx;
	double complex *vy;
	double complex *vz;
#ifdef BOUSSINESQ
	double complex *th;
#endif
#ifdef MHD
	double complex *bx;
	double complex *by;
	double complex *bz;
#endif
#ifdef COMPRESSIBLE
	double complex *d;
#endif
#ifdef WITH_PARTICLES
	struct Particle *part;
#endif
};

struct Grid {
    /* Total number of real grid points in each direction */
    int nxrtot;
    int nyrtot;
    int nzrtot;
    
    int nxtot;
    int nytot;
    int nztot;
    
    /* Number of cores in each direction */
    int nProcx;
    int nProcy;
    int nProc;
    
    /* current MPI rank in x and y*/
    int rankx;
    int ranky;
    
    /* Local number of complex grid points for current process in each direction */
    int nxloc;
    int nyloc;
    int nzloc;
    
    /* Local number of real grid points for current process in each direction
    CAUTION: only valid for non transposed transforms! */
    int nxrloc;
    int nyrloc;
    int nzrloc;
    
    /* start of the local array in the global array */
    int startx;
    int starty;
    int startz;
    
    int nrloc;  /* local number of real grid points */
    int nrtot;  /* total number of real grid points */
    
    int nloc;   /* local number of complex grid points */
    int ntot;   /* Total number of complex grid points */
                  
    
    double *kx, *ky, *kz, *kxt, *k2t, *ik2t;
    double kxmax, kymax, kzmax, kmax;
    
    double *mask;  
    
    int main;	/* Set to 1 only in one MPI process. Mark the main process (rank=0) */
};


struct VarName {
	int		length;		// Number of varnames
	char**	name;		// Varnames (need to be allocated properly)
};

struct Parameters {
	// Physics Parameters
	double lx;				/**< Box length in X*/
	double ly;				/**< Box length in Y*/
	double lz;				/**< Box length in Z*/
    
    int nx;                 /**< number of grid points in X*/
    int ny;                 /**< number of grid points in Y*/
    int nz;                 /**< number of grid points in Z*/
	
	double reynolds;		/**< Reynolds number (actully the inverse of the viscosity) */
	double reynolds_m;		/**< Magnetic Reynolds number (actully the inverse of the resistivity)  Used only when MHD is on*/
	
	double reynolds_th;		/**< Thermal Reynolds number (actully the inverse of the thermal diffusivity)  Used only when Boussinesq is on*/
	double reynolds_B;		/**< Reynolds number based on Braginskii viscosity */
	
	double x_hall;			/**< Hall parameter */
	
	double N2;				/**< Brunt Vaissala frequency squared */
	
	double omega_x;			/**< x component of the rotation vector. Only used when WITH_ROTATION is enabled. */
	double omega_y;         /**< y component of the rotation vector. Only used when WITH_ROTATION is enabled. */
    double omega_z;         /**< Vertical rotation rate (if Shear=1, Keplerian if found for (2.0/3.0). Only when WITH_ROTATION is on. */
    
	double shear;			/**< Shear rate (only when WITH_SHEAR is on) */
	
	double omega_shear;		/**< Pulsation of the time dependant shear (only when WITH_SHEAR and TIME_DEPENDANT_SHEAR is on, or alternatively WITH_LINEAR_TIDE) */
	
	double cs;				/**< Sound speed (only used when compressible is on) */
	
	// Particles parameters
	int    particles_n;		/**< Number of particles */
	double particles_mass;  /**< Mass of the particles */
	double particles_stime;	/**< Stopping time of the particles */
	double particles_dg_ratio;	/**< Dust to gas mass ratio for the particles feedback */
	double particles_epsilon;	/**< Pressure gradient epsilon */
	
	// Code parameters
	
	double cfl;				/**< CFL safety factor. Should be smaller than sqrt(3) for RK3 to be stable.*/
	double safety_source;	/**< Safety factor for SHEAR, Coriolis and Boussinesq terms (should be ~0.2 for accuracy) */
	
	double t_initial;		/**< Initial time of the simulation */
	double t_final;			/**< Simulation will stop if it reaches this time */
	double max_t_elapsed;	/**< Maximum elapsed time (in hours). Will stop after this elapsed time */
	
	int    interface_check;	/**< Number of loops between two checks for a user input. On slow filesystems, increase this number */
	int    interface_output_file;	/**< Set this option to force code outputs to a file instead of the screen */
	
	int    force_symmetries;	/**< set to enforce spectral symmetries and mean flow to zero. Useful when N^2 or kappa^2 < 0. (see enforce_symm() )*/
	int    symmetries_step;		/**< Number of loops between which the symmetries are enforced. Should be around ~20 for Boussinesq convection*/
	
	int    antialiasing;		/**< 2/3 Antialisaing rule. Could be removed if you assume is unimportant (untested feature). */
	
	int    restart;
	
	// Output parameters
	double toutput_time;		/**< Time between two outputs in the timevar file */
	double toutput_flow;		/**< Time between two snapshot outputs */
	double toutput_dump;		/**< Time between two restart dump outputs (restart dump are erased) */
	
	int		output_vorticity;	/**< Output the vorticity field in the 3D snapshots */
	
	struct VarName timevar_vars;	/**< Name of the variables needed in the timevar file */
	
	// initial conditions
	int	   init_vortex;			/**< Add a 2D Kida vortex in the box. Assumes S=1. Requires b>a*/
	double vortex_a;			/**< x dimension of the vortex */
	double vortex_b;			/**< y dimension of the vortex */
	
	int    init_spatial_structure;	/**< Init a user-defined spatial structure (see initflow.c) */
	
	int	   init_large_scale_noise;	/**< Init a large scale random noise */
	double per_amplitude_large;		/**< Amplitude of the large scale random noise */
	double noise_cut_length;		/**< Wavelength over which the noise is applied */
	
	int	   init_large_scale_2D_noise;	/**< Init a large scale 2D (x,y) random noise  */
	double per_amplitude_large_2D;		/**< Amplitude of the 2D large scale random noise */
	double noise_cut_length_2D;		    /**< Wavelength over which the 2D noise is applied */
	
	int    init_white_noise;		/**< Init a random white noise on all the modes */
	double per_amplitude_noise;		/**< total amplitude of the perturbation */
	
	int    init_mean_field;			/**< Force the mean magnetic field to a given value. */
	double bx0;						/**< Mean magnetic field in the x direction */
	double by0;						/**< Mean magnetic field in the y direction */
	double bz0;						/**< Mean magnetic field in the z direction */
	
	int    init_dump;				/**< Use a dump file as an initial condition (everything else, including t, noutput (...) is reinitialised) */
	
	int	   init_bench;				/**< Init the Benchmark initial conditions */
};

#endif

